/*************************************************************/
/*   Headers for a neo-hookean elasticity tensor in MOOSE    */
/*                    Status: compiles                       */
/*                                                           */
/* Copyright (C) 2018  eDgar                                 */
/*                                                           */
/* This program is free software: you can redistribute it    */
/* and/or modify it under the terms of the GNU General       */
/* Public License as published by the Free Software          */
/* Foundation, either version 3 of the License, or (at your  */
/* option) any later version.                                */
/*                                                           */
/* This program is distributed in the hope that it will be   */
/* useful, but WITHOUT ANY WARRANTY; without even the        */
/* implied warranty of MERCHANTABILITY or FITNESS FOR A      */
/* PARTICULAR PURPOSE.  See the GNU General Public License   */
/* for more details.                                         */
/*                                                           */
/* You should have received a copy of the GNU General Public */
/* License along with this program.  If not, see             */
/* <http://www.gnu.org/licenses/>.                           */
/*************************************************************/

#include "ComputeElasticityTensorNeohookean.h"
#include "MooseMesh.h"
#include "Assembly.h"

registerADMooseObject("spiderApp",
                      ComputeElasticityTensorNeohookean);

defineADValidParams
(
 ComputeElasticityTensorNeohookean,
 ADMaterial,
 params.addRequiredCoupledVar("displacements",
                              "The displacements appropriate"
                              " for the simulation geometry"
                              " and coordinate system");
 params.addParam<std::string>("base_name",
                              "Optional parameter that allows the user to define "
                              "multiple mechanics material systems on the same "
                              "block, i.e. for multiple phases");
 params.addClassDescription("Compute a neo-hookean"
                            "elasticity tensor in"
                            "Lagrange (referential)"
                            "coordinates."););

template <ComputeStage compute_stage>
ComputeElasticityTensorNeohookean<compute_stage>::ComputeElasticityTensorNeohookean(
    const InputParameters & parameters)
  :
  ADMaterial<compute_stage>(parameters),

  /* // Optional stuff (similar to */
  /* // ComputeElasticityTensorBase) */
  /* _base_name(// Test this condition */
  /*            isParamValid("base_name") */
  /*            ? */
  /*            // If test is true, return this */
  /*            adGetParam<std::string>("base_name") + "_" */
  /*            : */
  /*            // If test is false, return this */
  /*            ""), */
  /* _elasticity_tensor_name(_base_name + "elasticity_tensor"), */

  _ndisp(coupledComponents("displacements")),
  _grad_disp(3),
  /* _base_name(isParamValid("base_name") ? adGetParam<std::string>("base_name") + "_" : ""), */
  _elasticity_tensor(adDeclareADProperty<RankFourTensor>("elasticity_tensor"))
{
  // TODO:
  // if (adGetParam<bool>("use_displaced_mesh"))
  //   paramError("use_displaced_mesh", "The elasticity calculator needs to run on the undisplaced mesh.");

}

// TODO
// template <ComputeStage compute_stage>
// void
// ADComputeGreenLagrangeStrain<compute_stage>::
// computeQpProperties()
// {
//   issueGuarantee(_elasticity_tensor_name, Guarantee::ISOTROPIC);
// }

template <ComputeStage compute_stage>
void
ComputeElasticityTensorNeohookean<compute_stage>::
initialSetup()
{
  displacementIntegrityCheck();
  // fetch coupled variables and gradients (as stateful properties if necessary)
  for (unsigned int i = 0; i < _ndisp; ++i)
    _grad_disp[i] = &adCoupledGradient("displacements", i);

  // set unused dimensions to zero
  for (unsigned i = _ndisp; i < 3; ++i)
    _grad_disp[i] = &adGradZero();
}

template <ComputeStage compute_stage>
void
ComputeElasticityTensorNeohookean<compute_stage>::
displacementIntegrityCheck()
{
  // Checking for consistency between mesh size and length of the provided displacements vector
  if (_ndisp != _mesh.dimension())
    paramError(
               "displacements",
               "The number of variables supplied in 'displacements' must match the mesh dimension.");
}

template <ComputeStage compute_stage>
void
ComputeElasticityTensorNeohookean<compute_stage>::computeQpElasticityTensor()
{

  /* The elasticity tensor of the neo-Hookean model
   * in referential (a.k.a.) coordinates is
   *
   *  LaTeX:
   *  {\cal{C}}_{ijkl} = \lambda\,\C_{ij}^{-1}\,\C_{kl}^{-1}
   *                    +\mu\,\left(\C_{ik}^{-1}\,\C_{jl}^{-1}
   *                                +\C_{il}^{-1}*\C_{jk}^{-1}
   *                          \right)
   *
   *  - Belytschko, T., Liu, W. K., Moran, B., & Elkhodary,
   *    K. (2014). Nonlinear finite elements for continua and
   *    structures. : John Wiley and Sons, Ltd.
   *
   *  {\cal{C}}_ijkl = lambda * C_ij^-1 * C_kl^-1
   *                  + mu * (C_ik^-1 * C_jl^-1
   *                          + C_il^-1 * C_jk^-1);
   */

  ////////////////////////////
  // Mathematical constants //
  ////////////////////////////

  // You have to make sure that whatever operatations
  // (including addition, subtraction, etc. work as you
  // think). In this case, check that the definition of
  // RankFourTensor has the required operations.
  // [[moose/framework/src/utils/RankFourTensor.C]]
  /*
   * // LaTeX:
   * // \( \left(\delta_{ik}\,\delta_{jl}
   * //    + \delta_{il}\,\delta_{jk}\right) \)
   * //
   * // 1.0  0.0  0.0  0.0  0.0  0.0
   * // 0.0  1.0  0.0  0.0  0.0  0.0
   * // 0.0  0.0  1.0  0.0  0.0  0.0
   * // 0.0  0.0  0.0  0.5  0.0  0.0
   * // 0.0  0.0  0.0  0.0  0.5  0.0
   * // 0.0  0.0  0.0  0.0  0.0  0.5
   * // initIdentitySymmetricFour() in
   * // [[moose/framework/src/utils/RankFourTensor.C]]
   * //
   * // initIdentitySymmetricFour
   * // gives (de_ij de_kl + de_il de_jk)/2 -- andrew.wilkins
   * //
   * // It's possible that someone might use this with crazy
   * // non-symmetric stresses and strains and
   * // initIdentitySymmetricFour essentially symmetrises the
   * // strain before forming the stress. -- andrew.wilkins
   * //
   * // See also [BelytschkoLiuMoranElkhodary14, p.240]
   * const RankFourTensor dik_djl__dil_dkj = RankFourTensor::initIdentitySymmetricFour;
   */
  // Identity matrix
  // 1  0  0
  // 0  1  0
  // 0  0  1
  // initIdentity() in
  // [[moose/framework/src/utils/RankTwoTensor.C]]
  const RankTwoTensor dij(RankTwoTensor::initIdentity);

  /*
   * // LaTeX:
   * // \( \delta_{ij} \, \delta_{kl} \)
   * //
   * // 1  1  1  0  0  0
   * // 1  1  1  0  0  0
   * // 1  1  1  0  0  0
   * // 0  0  0  0  0  0
   * // 0  0  0  0  0  0
   * // 0  0  0  0  0  0
   * const RankFourTensor dij_dkl = dij.outerProduct(dij);
   *
   * // // Another way for dij_dkl
   * // // To make sure that nothing is out of place, fill the
   * // // whole tensor with zeros
   * // // [[moose/framework/include/utils/RankFourTensor.h]]
   * // dij_dkl.zero();
   * // // Set the corresponding values
   * // // [[moose/framework/src/utils/RankFourTensor.C]]
   * // dij_dkl.fillGeneralIsotropic(1, 0, 0);
   */

  /////////////////////////
  // Material parameters //
  /////////////////////////

  // Hard-code Young's modulus for testing
  const Real E = 30e6;
  // Hard-code Poisson's ratio for testing
  const Real nu = 0.4;

  // Value of Lamé constants in terms of Young's modulus
  // and Poisson's ratio
  const Real mu = E / (2.0 * (nu + 1.0));
  // // 3 multiplications, 1 division, 1 subtraction, less clear:
  // _lambda = 2.0 * _mu * _nu / (1.0 - 2.0 * _nu);
  // 3 multiplications, 1 division, 1 addition, 1 subtraction:
  const Real lambda = E * nu / ((nu + 1) * (1 - 2 * nu));

  /////////////////////////
  // Kinematic variables //
  /////////////////////////

  for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
  {
    // Displacement gradient
    ADRankTwoTensor dudX((*_grad_disp[0])[_qp],
                         (*_grad_disp[1])[_qp],
                         (*_grad_disp[2])[_qp]);
    // Deformation gradient tensor F = 1 + du_{i}/dX_{j}
    // Just copy the value for the given quadrature point
    ADRankTwoTensor F = dudX + dij;

    // Inverse of the Right Cauchy-Green deformation tensor
    ADRankTwoTensor Cij_1 = (F.transpose() * F).inverse();

    /* // Why doesn't this work? */
    /* ADRankFourTensor C_ijkl = lambda * Cij_1.outerProduct(Cij_1) */
    /*   + (mu - lambda) * (Cij_1.mixedProductIkJl(Cij_1) */
    /*                      + Cij_1.mixedProductIlJk(Cij_1)); */

    _elasticity_tensor[_qp] = lambda * Cij_1.outerProduct(Cij_1)
      + (mu - lambda) * (Cij_1.mixedProductIkJl(Cij_1)
                         + Cij_1.mixedProductIlJk(Cij_1));
  }
}
